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ABSTRACT 

In this paper we study the shock/longitudinal vortex interaction problem in axisymmet- 
ric geometry. Linearized analysis for small vortex strength is performed, and compared with 
results from a high order axisymmetric shock-fitted Euler solution obtained for this purpose. 
It is confirmed that for weak vortices, predictions from linear theory agree well with results 
from nonlinear numerical simulations at the shock location. To handle very strong longitu- 
dinal vortices, which may ultimately break the shock, we use an axisymmetric high order 
essentially non-oscillatory (ENO) shock capturing scheme. Comparisons of shock-captured 
and shock-fitted results are performed in their regions of common validity. We also study 
the vortex breakdown as a function of Mach number ranging from 1.3 to 10, thus extending 
the range of existing results. For vortex strengths above a critical value, a triple point forms 
on the shock, leading to a Mach disk. This leads to a strong recirculating region downstream 
of the shock and a secondary shock forms to provide the necessary deceleration so that the 
fluid velocity can adjust to downstream conditions at the shock. 


'Research supported by NASA contract NAS1-19480 while the authors were in residence at ICASE, 
NASA Langley Research Center, Hampton, VA 23681-0001. 

2 ICASE. Mail Stop 132C, NASA Langley Research Center, Hampton, VA 23681-0001. 

3 Program in Computational Science and Engineering, Florida State University, Tallahasee, FL 32306 
4 Division of Applied Mathematics. Brown University, Providence, RI 02912. Research of this author was 
also supported by ARO grant DAAH04-94-G-0205, NSF grant DMS-9500814, NASA grant NAG1-1145 and 
AFOSR. grant 95-1-0074. 


l 




1 Introduction 


Over the last 15 years, there have been numerous experimental (Cattafesta k Settles 1992; 
Cattafesta, 1992; Delery, Horowitz. Leuchter k Solignac 1984; Dosanjh k Weeks 1965; Nau- 
mann k Hermans 1973), theoretical (Chang 1957; Ribner 1954), and computational (Ellzey, 
Henneke, Picone k Oran 1995; Erlebacher, Jackson k Hussaini 1996; Kopriva 1988; Las- 
seigne, Jackson k Hussaini 1991; Meadows 1991; Meadows, Kumar k Hussaini 1991; Mead- 
ows 1995; Pao k Salas 1981; Rizetta 1995; Zang, Hussaini k Bushnell 1984) studies of the 
shock-vortex interaction problems. Major effort has been devoted to investigating the inter- 
action of a shock with either plane waves (Ribner 1954; Ribner 1986; McKenzie k Westphal 
1968; Chang 1957; Erlebacher et al. 1996; Lasseigne et al. 1991) or with a vortex w r hose axis is 
aligned with the shock (Naumann k Hermans 1973; Ellzey et al. 1995; Erlebacher et al. 1996; 
Chang 1957; Pao k Salas 1981; Zang et al. 1984). In these studies, the main interest was 
to understand the shock distortion, the vorticity amplification mechanisms engendered by 
the shock, and the properties of the radiated noise. The observations of the experiments 
(Naumann k Hermans 1973; Cattafesta 1992; Delery et al. 1984; Dosanjh k Weeks 1965), 
have been in general substantiated by the theoretical work of Ribner (1954) and Ribner 
(1986), based on linear theory. Numerical simulations have recently begun to quantify some 
of the experimental results, in both the linear and nonlinear regimes (Ellzey et al. 1995; 
Erlebacher et al. 1996; Meadows et al. 1991: Meadows 1995). 

The configuration in which the vortex has its axis normal to the shock occurs in practice, 
for example, when a wing tip vortex shed from a canard of a supersonic fighter plane intersects 
the shock that lies over the wing. The resulting deceleration of the vortex can lead to vortex 
breakdown (if the vortex and shock strengths are appropriate), resulting in destabilizing 
forces on the airplane (Delery et al. 1984). Thus, a theoretical, or heuristic determination of 
the conditions under which this breakdowm occurs is of practical interest for both the design 
of fighter planes, and for the problem of controlling the resultant destabilizing forces. 

There have been only a handful of experiments concerned with the longitudinal vor- 
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tex/shock interaction (Cattafesta k Settles 1 992 5 Cattafesta 1992; Delery et al. 1984). In 
these cases, particular care was taken to ensure that the vortex was axisymmetric, and that 
it interacted with a planar normal shock. The experiments differed in the shock creation pro- 
cess. For example, Cattafesta k Settles (1992) and Delery et al. (1984) created a shock using 
a 2-D pitot type air intake mounted in a section of uniform flow. In addition, Cattafesta 
k Settles (1992) also created a normal shock wave by over-expanding the exit flow from a 
supersonic nozzle. The associated numerical simulations were able to reproduce the gross 
features of the experimental results. 

To enhance our understanding of the shock-induced vortex breakdown phenomena, we 
study the case of a shock of infinite extent interacting with a longitudinal vortex. We find 
that the flow is not always steady. However, for moderately weak vortices, a steady state 
appears possible. The vortex breakdown is known to be a function of the helix angle which 
is defined as the arctangent of the ratio of the maximum azimuthal velocity and the mean 
axial velocity. In the incompressible case, it is observed that once the helix angle attains 
the critical value of 57°, the vortex is prone to breakdown (Delery et al. 1984; Spall, Gatski 
k Grosch 1987). The critical helix angle is lower for compressible flows. In the presence 
of a shock, the critical value of the helix angle decreases further due to the higher flow 
deceleration. It is attained more easily for stronger vortices and higher shock Mach numbers 
(Cattafesta 1992; Delery et al. 1984). If the vortex circulation is further increased, strong 
nonlinear effects come into play. The pressure associated with the vortex core, which scales 
quadratically with circulation, leads to nonlinear effects responsible for the formation of a 
Mach disk along with a strong downstream recirculating zone with a complex structure. As 
in the previous work by Erlebacher et al. (1996) where the vortex axis was parallel to the 
shock, nonlinear effects are directly related to the product of vortex strength and shock Mach 
number. This is the direct result of a linear scaling of vortex strength with circulation, and 
a quadratic scaling of the maximum pressure variation within the vortex core. 

The paper is organized as follows: Section 2 contains the mathematical formulation of the 
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problem. Section 3 provides details of the linear analysis, including a simplified high Mach 
number expansion of the Rankine-Hugoniot conditions. Section 4 describes both the shock- 
fitted compact (S-F) and the shock-capturing essentially non-oscillatory (ENO) numerical 
algorithms. Section 5 presents the consistency checks of both algorithms against linear 
theory and against each other. In Sections 6 and 7, we study the influence of the variation 
in vortex strength and shock Mach number on vortex breakdown and shock bifurcation. We 
first compute the vortex breakdown curve as a function of shock Mach number ranging from 
1.3 to 10. Then, in Section 7, we consider vortex strengths leading to the formation of a 
Mach disk and strong recirculating zones downstream of the shock. Unsteady effects and the 
influence of initial conditions are also addressed in this section. Some concluding remarks 
are made in the final section. 

2 Problem Formulation 

We seek to study the interaction of a longitudinal vortex with an infinite shock. For sim- 
plicity, the vortex is assumed to be axisymmetric with its axis perpendicular to the shock of 
infinite radial extent. The geometry is shown in Figure 1. 


Freestream 



Figure 1: Geometry. 
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The flow is assumed to be governed by the axisymmetric compressible Euler equations 
(in conservative form): 


dp d{pu x ) d[pUr) pur 

dt dx dr r 

djpUx) d{pu 2 x +p) d(/9U r U r ) pUxUr 

dt dx dr r 

djpUr) djpUrll x ) d{pu 2 r + p) p{u l - tig) 

dt dx dr r 

d{pue) d{pueu x ) d{pugu r ) 2pugu r 

~~dt~ + dx + d~r + r 

dE d(u x (E + p)) d(u r (E + p)) + p) 

~di + dx + dr r 


( 2 . 1 ) 


where p is the density, (ii x .u r .u$) is the velocity vector in the axial, radial, and azimuthal 
directions, and E is the total energy, 


E = — —r + \p{u 2 x + u 2 r + u l) ( 2 - 2 ) 

7 — 1 2 

with the ratio of specific heats 7 = 1.4 for air. The pressure, density, and temperature are 
nondimensionalized with respect to their mean upstream values Pi, pi, and T\ respectively 
and are related by the ideal gas law 

p = pT. (2.3) 


The velocity is scaled by the reference velocity c* = \fT[, related to the upstream mean 
sound speed Ci = . In the frame of reference in which the mean shock is stationary, the 

mean upstream Mach number is therefore M x = \U X \/ y/j, where M x is the upstream Mach 
number in a frame of reference in which the mean shock is stationary. Note that the chosen 
nondimensionalization leaves the Euler equations (2.1) invariant. Finally, lengths are scaled 
by the vortex core radius r 0 . 

The shock is initially located at x = 0 plane; the axial extent of the upstream domain 
extends to x = b. while the leftmost boundary on the downstream side is located at x = a. 
The flow is from right to left so that a < 0, and b > 0. The axis of symmetry is r = 0, and the 
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radial coordinate extends to r = c, where c is sufficiently large so that all flow disturbances 
are effectively zero at this freestream boundary. 

At t = 0, we choose a mean flow consistent with a stationary shock at x — 0. Upstream 
of the shock (*>0), 


P= 1, U\ = Ux = —\/lM, U r = Ug = 0, P, = p = 1 , (2.4) 


_ ( 7 + 1 W? 

r ~ (7— l)Jl/j i +2 ' 

U r = Ug = 0. 


(2.5) 


while the downstream mean solution is (.r < 0), 

u 7 = u = 

t2 x ( 7+1 )Mi ’ 

P, = n = glMnlllzll 

Next, we superimpose on the mean flow an isentropic vortex with its axis along r = 0. 
Analytical forms of such vortices which are steady-state solutions of the Euler equations 
exist with arbitrary radial profiles. We choose an exponentially decaying profile to reduce 
interactions of the vortex with the outer domain boundary in the radial direction. The 
perturbations of azimuthal velocity u f $ and temperature T* due to the vortex are given by 


u e 


■ = H e 0.S(.-r’) T= J_ 


( 2 . 6 ) 


2tt 877r 2 ro 

where r 0 is the vortex core radius and t is a nondimensional circulation at r = 1, related to 
the dimensional circulation T by 


r 0 c 


(2.7) 


The axial and radial velocities u ' x , u' r and the perturbation entropy S' = log {pjp 1 ) are 
constant inside the vortex; Ug is maximum at r = 1. Because of the particular radial profile 
of u'g , the vortex circulation decays exponentially fast to zero as r — » oo. Another measure 
of the relative strength of the vortex is given by the ratio of the maximum u' e to U\. This 
quantity, denoted by r, measures the inclination of the streamlines with respect to the 
symmetry axis of the vortex (Delery et al. 1984). It is also the tangent of the helix angle. 
The relationship between e and r is 


e 

27tM v /7 


( 2 . 8 ) 
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The upstream vortex defined by (2.6) has the following features: 

• It is a steady-state solution of the Euler equation (2.1) if u x = 0. Otherwise, it is 
passively convected by the flow with velocity u x , 

• it has constant entropy, 

• it is regular at the symmetry axis, 

• it decays exponentially fast as r -> oo, thus ensuring that freestream boundary effects 
are virtually absent. In fact, in most calculations we set the computational boundary 
at r = 5, at which the relative effect of the perturbation on the pressure and density 
is on the order of 10 -l °. 

Previous studies (Cattafesta 1992; Delery et al. 1984) have indicated that for low and 
moderate vortex strengths, a steady state is eventually reached. However, the shocks in 
these simulations were kept stationary by either placing the shock in a divergent nozzle 
(Delery et al. 1984), or by attaching a Mach disk at the end of an underexpanded jet 
(Cattafesta 1992). Such artifacts accelerate convergence to a steady state. Should a steady 
state exist, one often assumes that the final state is independent of initial conditions. For 
example, Delery et al. (1984) assume that the total enthalpy is independent of time. They 
therefore simplify their calculations by building this constraint into their equations. However, 
this precludes even genuinely unsteady solutions. In contrast, we do not presume that steady 
solutions exist, but compute them as they evolve in time. 

To lessen possible transients during the initial stages of the simulation, we multiply the 
initial conditions by a function s(;r). We have considered three different functions: 


IC-a 


s(x) = 1.0 


(2.9) 


IC-b 


s(z) = 


( 


0 

1 


x < 0 
x > 0 


( 2 . 10 ) 
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IC-c 


x < a 


s(.r) = 



a < x < f3 
x > 0. 


( 2 . 11 ) 


In case IC-a, the vortex structure (defined by u'e and T') is independent of the axial location 
and exists across both upstream and downstream domains. Then a shock is externally 
generated along x = 0 at t = 0. The initial conditions satisfy steady-state Euler equations, 
but the pressure and density jumps across the initial shock do not satisfy the steady Rankine- 
Hugoniot conditions, (which then causes the shock to adjust itself.) In case IC-b, the vortex 
is only defined in the upstream domain, and abruptly ends at the shock, and it offers a clean 
model configuration. Finally, case IC-c is intermediate between the first two cases. The 
transition function smoothly varies from zero to one and is only nonzero in the upstream 
domain. Thus, the Rankine-Hugoniot conditions are initially satisfied, but the upstream 
modified vortex is no longer a solution to the steady-state Euler equations. Note that IC-b 
is a special case of IC-c where a — > 0, 0 — ¥ 0. Numerical tests presented in a later section 
will show that IC-b and IC-c actually produce similar transient structures, albeit shifted in 
time. 

The boundary conditions are chosen as follows. We prescribe all variables at the super- 
sonic inflow on the upstream side (a: = b). At the subsonic downstream boundary, we employ 
the characteristic boundary conditions for the ENO method, or a buffer domain technique 
for the S-F method. The details are discussed in Section 3. The freestream boundary (r = c) 
is sufficiently removed from the vortex so that a simple Neumann boundary condition is suf- 
ficient for the ENO algorithm. The shock-fitted (S-F) algorithm permits a choice between 
characteristic conditions and Dirichlet boundary conditions. Naturally, w r e impose symme- 
try conditions at r = 0. Radial derivatives of all variables are set to zero, except for the 
azimuthal and radial velocities which vanish. 



3 Numerical Methods 


We adopt two numerical algorithms of high order accuracy to solve the unsteady axisymmet- 
ric compressible Euler equations for the shock/vortex interaction problem — 1 ) a shock-fitted 
(S-F) algorithm and ii) an essentially nonoscillatorv (ENO) algorithm which captures the 
shock. The former is based on a sixth order compact scheme and has the advantage of 
precluding Gibb’s phenomena associated with shock capturing; however, its disadvantage 
is that it cannot handle very weak shocks or situations of strong shocks with triple points. 
The latter is based on adaptive stencil interpolation and nonlinearly stable time discretiza- 
tion. which ensures formally third order accuracy in the smooth part of the solution while 
maintaining sharp, essentially nonoscillatory shock transitions. 


3.1 Shock-fitting method 


We provide here a brief description of the numerical method (see Erlebacher et al. (1996) 
for more details). All spatial derivatives in the shock-fitting method are computed with 
a compact 6 th order discretization scheme (Carpenter, Gottlieb & Abarbanel 1993) which 
is stable for Dirichlet boundary conditions. Both x and r coordinate directions are non- 
periodic. Derivatives at the first and second points at both ends of the computational 
domain are explicitly defined by 

7 

hu 0 = ^ d{ U{ 

i = 0 
7 

hlIXy — 6, u, 

1=0 

where the parameters a, and b x are given by 


296 415 125 985 215 791 25 145 

^ 105 ’ Is"’ IT’ ~ 48 ' ~l 2 ' 336 * 

3 _ 211 109 115 _1 _ 23 _ _J_ 

* l 6 ’ _ 180 ’ 48 24 ’ 144 ’ 3 ’ 240 ’ 72 ^’ 


( 3 . 1 ) 

( 3 . 2 ) 


and h is the uniform grid spacing. Identical explicit formulas are derived for u' N and u' N _ 1 , 
where N is the number of grid points along the nonperiodic direction. This stencil is com- 
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monly referred to as 5 2 — 6 — 5 2 (Carpenter et al. 1993). The derivatives at interior nodes 
(z = 1, • • • , N — 1) are computed via the standard scheme (Lele 1992) 

J2 = /r ‘ Pj u >+j (3-3) 

.? = -! j =- 2 

with coefficients q±i = 1. Oo = 3, 0±2 — 1/12, /?±i = 7/3 and 0 O — 0. 

To maintain the Q th order accuracy of derivatives in the radial direction at the axis, we 
build in the required symmetry properties into the matrix derivative operator. This requires 
modifying the formulas for the first two radial points which are now solved as part of the 
implicit system. Let Uj be the value of any variable at the j th radial point, with j = 0 
corresponding to the symmetry axis. Symmetry conditions imply that r<_j = uj, while 
antisymmetric conditions imply u_j = — uj. Combining these requirements with the interior 
formula (3.3), we obtain for the point on the axis under symmetric conditions. c*o = 3, 
Qj = 0, 3 t = 0. The first point off the axis has 02 = ~0o = 7/3, 03 = —0\ - 1/12, 
while the a's take on the interior values for the first point off the axis. The requirement of 


antisymmetry leads to Of 0 — 3, Qi = 2, 0q = 0, 0\ — 14/3, 02 = 1/6, 03 = 0 for the axis 
point, while 0\ = —0-\ — 7/3, 0o = 02 — 1/12, and once again, the a's take on the interior 
values, for the first point off the axis. 


The symmetry axis is a geometric singularity, and requires special treatment. We rec- 


ognize that it is an interior point, and the flow variables satisfy the Euler equations on the 
axis. We therefore impose the Euler equations, where singular terms have been evaluated 


using L’Hospitale rule. For example, u r /r is replaced by du r jdr at the axis. 
The grid is stretched in the radial direction according to 

"sinh (V — Yq)0 


r = r 0 


sinh 0Yo 


+ 1 


(3.4) 


w’here 


1 . 1 + (e* - l)(r 0 /r max ) 

2 A ° g 1 + (e~ A - 1 ){r 0 /r max ) 


(3.5) 


w-here 0 < V < 1 is the computational coordinate in the radial direction. 


This choice of 


stretching allows the grid points to concentrate around r = 1 and extend to a maximum 
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radius of r max ; the stretching of the grid is controlled by the parameter A. As detailed in 
Erlebacher et al. (1986) we transform the downstream physical axial coordinate x using the 
coordinate transformation 

A' = : T ~ a - (3.6) 

x s (r,t) - a 

where :r s (r, t) is the shock shape, and x = a is the leftmost boundary of the downstream 
domain. 

At the downstream boundary, we implement the buffer domain technique introduced 
by Ta’assan k Nark (1995). We modify the downstream Euler equations in a finite buffer 
domain abutting the downstream boundary, by replacing d/dt by d/dt + Uq(X)V, to control 
the direction of the characteristics and the characteristic speeds. At the exit plane, X = 0. 
the characteristic speeds become 112 (b) + Uo(b), Uo(b) + l'o(b) ± C 2 (&). The axial profile for 
U 0 (x) is a power function of x described in Erlebacher et al. (1996). 

We briefly present the derivation of a time evolution equation for the shock motion 
(Hussaini, Kopriva, Salas k Zang 1984, Canuto, Hussaini, Quarteroni k Zang 1987; Kopriva 
1991). Given the shock shape, x = x s (r, t) and the Euler equations in the form u f = R u and 
p t — R the characteristic normal to the shock pointing in the upstream direction has the 


( P.t - R p ) + (n- (u t - R u ) = 0 (3.7) 

where £ = P 2 V 1 T 2 is a function of the downstream temperature and density, and n is the 
unit normal to the shock pointing in the upstream direction: 

n= yrrif:) • (3 - 8 » 

R u and R p are numerically computed as part of the Runga-Kutta advancement. To apply 
Eq. (3.7) at the shock requires the time derivatives of pressure and velocity. These are ob- 
tained by solving the Rankine-Hugoniot conditions the downstream pressure p 2 and velocity 


2 , 7-1 f 2 \ 

P 2 = Pi+M) 
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(3.9) 


X 7 - 1 r , 27 Pi 

■ 2 — — -r*i-( — 7 — Y' 

7 + 1 7 + 1 pi^i 

and differentiating {p 2 ,S 2 ) with respect to time. Here. Si = u,n — u s , where u s = w s x-n is the 
shock velocity normal to the shock, while i and w s are the unit vector and the instantaneous 
shock velocity respectively. The required time derivatives are readily computed and inserted 
into the characteristic equation above, to obtain the final evolution equation 


(n x iv s )' t 


( F - C-4)(n • uj),* + C(u 2 • n |( ) + (D - (B)p ut + (E - SC)p u - H 

F + C(l-A) 


(3.10) 


w r here 

'F — R-p (, n 2 • R u , 


(3.11) 


and the constants .4, F . D, B , E : C are functions of the upstream flow variables, and the 
shock velocity (Erlebacher et al. 1996). 

Time advancement is based on a low-storage, 5 stage, 4 th order accurate Runga-Kutta 
scheme (Carpenter & Kennedy 1994). Although 5/3 more expensive than Williamson's 3rd 
order low-storage scheme at fixed CFL, the achievable CFL is higher by a factor of 1.9, w r hich 
is greater than 5/3, so with no change in memory requirement, there is a net gain of CPU 
time, along with the increased accuracy. The interior equations and shock motion equations 
are advanced simultaneously in time. The Rankine-Hugoniot conditions are applied at the 
end of each Runga-Kutta stage to update the downstream variables at the shock. 


3.2 Shock-capturing method 

The shock capturing method w'e use in this work belongs to the class of high order nonlinearly 
stable essentially nonoscillatory (ENO) methods developed by Shu k Osher (1988): Shu k 
Osher (1989); and Shu. Zang, Erlebacher, Whitaker k Osher (1992). The detailed description 
of the algorithm, along with information on efficient implementation can be found in these 
references. Here we only highlight a few key points and describe issues which relate to the 
applications of ENO schemes to the particular system (2.1). 
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The main idea underlying an ENO scheme lies at the approximation (interpolation) 
level. For piecewise smooth functions and a fixed stencil, high order interpolation schemes 
inevitably cross the discontinuities, causing not only loss of accuracy but also over- and under- 
shoots (Gibbs phenomenon). ENO (Harten, Engquist, Osher Sz Chakravarthy 1987) is an 
adaptive stencil interpolation, with the local stencil chosen as the smoothest possible among 
all candidates. This, together with upwinding (realized through flux splitting), characteristic 
decomposition (which effectively decouples the system locally), and a nonlinearly stable time 
discretization, ensures that the scheme achieves high order resolution in the smooth part 
of the solution while maintaining sharp, non-oscillatory shock transitions. For the shock- 
vortex interaction problem considered in this paper, which contains both strong shocks and 
complex structures in the smooth part of the solution, a high order ENO scheme is an ideal 
candidate. We use the uniformly third order ENO scheme (fourth order in Li), based on 
a Lax-Friedrichs building block, and a third order total-variation-diminishing Runge-Kutta 
time discretization (Shu et al. 1992). 

The procedure underlying ENO schemes is best described with reference to the one- 
dimensional scalar version of (2.1): 


du dfju) _ 
dt dx 


(3.12) 


and then indicate how it can be generalized to solve the full system (2.1). The spatial 
derivative in (3.12) is discretized by the conservative difference 


dfju) 

dx 


Ax 


(/i+i-A-i)- 


A linear combination of the point values of the flux function /(«), 

fcl+r+l 

fj+ 1= E C ( k h k )f( U J+k) 

k=k] 


(3.13) 


(3.14) 


is used to represent the numerical flux. Here, r is the order of accuracy in L°° (i.e., the 
scheme is (r -f l) (/l order L l ), c is a constant matrix independent of /(u), hence is computed 
only once and stored (precise formulas are found in Shu et al. (1992), and the leftmost point 
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location k\ is decided locally by the ENO interpolation procedure (the exact procedure is 
detailed in, e.g., (Shu et al. 1992), and can differ from point to point.) Upwinding is used 
in this stencil choosing process. 

With an explict ENO scheme, a multi-dimensional problem scalar equation is handled 
dimension by dimension. Each derivative is treated as one-dimensional in computational 
space, and the above algorithm is applied. For systems of multidimensional equations, a local 
characteristic decomposition is first performed to transform the equations into a decoupled 
set of multidimensional scalar equations whose unknowns are the Rieman invariants of the 
original system of equations (Shu et al. 1992). Derivatives of the Rieman invariants are 
computed according to the one-dimensional algorithm. 

The ENO scheme treats the symmetry axis differently from the compact scheme algo- 
rithm. The point value ENO procedure stores solution variables at the cell centers and 
computes fluxes at cell edges. The axis r = 0 is placed at a cell edge rather than at a cell 
center. Therefore, the computational domain is extended to the other side of the axis and 
uses of ghost points. Symmetry conditions are imposed on all the variables at the end of 
each iteration at the ghost points, except for the radial and tangential velocities which are 
antisymmetric about r = 0. 

4 Linear Analysis 

This section examines the linearized Euler and Rankine-Hugoniot (R-H) conditions with 
a view to obtaining some analytical and physical insight into the behaviour of the flow. 
Consider a shock normal to the mean flow and label the upstream and downstream sides by 
subscripts 1 and 2 respectively. From the steady-state Rankine-Hugoniot (R-H) conditions, 
one readily deduces the relations 

_ 2 + Af?( 7 -P 

2 7 M 1 2 -(7-l) 

= (7 + 1 )Mj 

P 21 Pi 2 + Mj 2 ( 7 — 1) 
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(4.1) 


P± = 2 jMf - (7 ~ 1 ) 

Pi 7 + 1 

between mean upstream and downstream variables. 

Now consider an infinitesimal perturbation wave (denoted by primes) superimposed on 
the mean flow, but not necessarily aligned with it. These perturbation velocities are decom- 
posed normal and parallel to the mean shock position. This in turn induces a perturbation 
in the shock position and orientation. Together with an assumed shock shape given by 

x = x s {r,t), (4.2) 

and the introduction of a finite shock velocity into the R-H conditions, their linearization 
with respect to perturbation amplitude leads to a linear system of equations for the down- 
stream variables. For simplicity, we set the upstream axial and radial perturbation velocity 
components, and the upstream perturbation entropy, to zero. Solving the linearized R-H 
conditions for the downstream perturbation variables yields 


£2 

= kl2 — 

„ X s t 

+ n x -^ 

(4-3) 

7 

7Pi 

Cl 


P 2 

\ p i 
= A 22 

„ x st 

+ n 2 — 

(4.4) 

IP2 

IPi 

Cl 



_ A Pl 
— A 32 

+ n 3 — 

(4.5) 

c 2 

7Pi 

Cl 


K 2 

— H 4 X Sfr . 


(4.6) 

C 2 





Here — * is the downstream perturbation entropy. 

1 P2 P2 

The nomenclature for the matrix coefficients is taken from Chang (1957) where the gen- 
eral solution for arbitrary upstream disturbances was derived. With the previous definitions, 
the matrix coefficients are 

a . 2 = + 'frjjj2 {(i -fti)! 1 + (i'- OMfl + [] - (^) 2 )} <0 
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A 22 = — ~ ~ -Jp ~ ^2l) + (1 - ^2 H 1 + (7 ~ 1)(1 “ ^2l)M|]| > 0 

yV/2 

A32 = ,v/ ( 1 ^ M 2 ) (1 — P21 ){2 + (7 _ 1 )( 1 _ A21 ) Af 2 2 } < 0 

n, = (7 - i) (1 ~-^ 21 £ a/ 2 > 0 

P21 

n 2 = — [2 + (7-l)(l-«i)Af|]>« 

1 - A/ 7 P21 

n.3 = -1-^72— ~ [* + ^1 + (7 -!)(!- ^2i)M 2 2 ] >0 

1 — M 2 P21 

n 4 = -M 2 (P21-1) <0. (4.7) 

The inequalities in the expressions above hold for the range of upstream Mach numbers 
Mi > 1 Several conclusions can be drawn from Equations (4.3). First, note that u' g is con- 
tinuous across the shock, and is therefore absent from the linearized R-H relations. Second. 
u' r2 has the same radial profile as the shock slope z s , r (r, t). If the shock moves in the positive 
x direction, then x sj > 0. In addition, if the axial velocity decreases aw'ay from the axis, 
then x s _ r < 0, which implies that u' r2 > 0. The sign of u' x2 is also positive since p\ < 0 in the 
upstream vortex core. Thus, if the shock is displaced in the upstream direction, its velocity 
is relatively fast at the axis, and a counter-clockwise (in region r > 0) downstream vortex 
ring is generated. This vortex ring then convects downstream with the mean axial velocity 
t'2, while new vortex rings are continuously generated by the steady upstream. The final 
result is a perturbation velocity u' x2 pointing upstream at the axis, and pointing downstream 
at the top of the vortex ring. Superimposed on this ring of azimuthal vorticity is a down- 
stream axial vortex characterized by the same azimuthal velocity distribution as upstream 
(to leading order), but with a lower relative (with respect to the downstream mean pressure) 
pressure variation within the vortex (0.67 < A 2 2 < 0.72, for 1 < M < oc). As expected, the 
signs of A12 and II 1 are consistent with positive entropy generation across the shock. 

Several numerical simulations of shock/vortex interactions (with the vortex axis parallel 
to the shock), have shown that a wave propagates along the shock at a velocity different 
from either the upstream or the downstream sound speed (Erlebacher et al. 1996; Meadows 
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Figure 2: Propagation speed of disturbances along the shock in the linear regime, normalized by 
downstream sound speed. 

1995). The acoustic wave, created by the shock perturbation, expands outward from its 
source at the downstream mean sound speed. The center of the instantaneous acoustic pulse 
convects downstream at velocity U 2 . Its speed along the shock is easily calculated to be 
C* = C 2 yfl - A/f- 

5 Discussion of Physical Results 

In this section, results from the S-F and ENO algorithms, as applied to the problem of 
shock/vortex interaction, are compared to the predictions of linear theory and against each 
other. Because of the low-order accuracy of ENO results in the vicinity of the shock, detailed 
verification of shock-related information is rather difficult. However, it is possible to verify 
results of ENO scheme in the downstream region away from the shock where pointwise or 
integrated quantities can be checked. Therefore, after providing a mutual verification of 
S-F method and linear theory in the linear regime, we use the S-F results to establish the 
limitations of linear theory. Then we use the S-F results to validate the ENO results for 
moderately strong vortices, especially with reference to the structure of the flow away from 
the shock. Finally we venture with the ENO scheme into the highly nonlinear regime of 
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stronger vortices, which cannot be handled by the S-F method. 

5.1 Linear regime 

The validity and accuracy of the S-F algorithm have been verified in numerous situa- 
tions involving plane acoustic, shear or entropy waves (Zang, Kopriva & Hussaini 1983, 
Zang et al. 1984, Lasseigne et al. 1991, Erlebacher et al. 1996). We consider specifically ax- 
ial vortices (2.6) with strengths in the range of e € [10' 3 ,3] interacting with a M = 2 shock. 
After the downstream solution computed by the S-F method has evolved for a fixed period 
of time, the downstream values of field and thermodymamic variables are extracted at the 
shock. These results are then compared with the corresponding predicted values obtained 
from the linear equations (4.3)-(4.6) using the upstream values and the shock speed at the 
particular instance in time. 

We compute the error (both absolute (a.e.) and relative (r.e.)) for entropy, pressure and 
velocity normal to the shock. Results at the symmetry axis (where the maximum error is 
found to occur) are presented in Table 1 and 2. As expected, the linear theory prediction 
is found to be accurate and match the shock-fitted results in the range of vortex strength 
e < 1 (r < 0.067). 

Since the perturbation p[ = <3(e 2 ), u\ x = 0, and u' e is continuous across the shock, it 
follow's that all downstream perturbation variables are 0{e 2 ). The absence of 0(e) terms 
implies that the absolute error is 0(e 4 ) in the linear regime, and the relative error is 0(e 2 ). 
These results are clearly confirmed in Table 1, except for where absolute errors reach the 
roundoff error of the computer (< 10~ 13 ). 

A threshold of nonlinearity was previously established in Erlebacher et al. (1996) where 
it was deduced that nonlinear effects become significant when Me > 2, at least when the 
vortex axis is parallel to the shock. Although this derivation does not strictly hold when the 
vortex axis is perpendicular to the shock, Table 2 does indicate that relative errors related 
to nonlinear effects grow beyond 1% when e is in the 1 to 2 range. Note that by e = 3 
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t 

i 

0.1 

0.01 

0.001 

T 

0.067 

0.0067 

0.00067 

0.000067 

S (a.e.) 

5.8 x 10 -5 

5.8 x 10~ 9 

9.1 x 10” 13 

2.1 x lO" 14 

S (r.e.) 

2.4 x 10"* 

2.4 x 10' 4 

3.7 x 10" e 

8.5 x 10“ 6 

p (a.e.) 

5.3 x 10~ 4 

5.1 x 10“ 8 

2.7 x 10" 12 

7.8 x 10~ 14 

P ( r - e -) 

4.9 x 10" 3 

4.6 x 10" 5 

2.5 x 10 -t 

7.0 x 10“' 

u n (&.e.) 

6.0 x 10~ 5 

6.1 x 10~ 9 

9.1 x 10" 13 

3.1 x 10~ 15 

(r.e.) 

1.3 x 10" 2 

1.4 x 10" 4 

2.1 x lO’ 6 

7.2 x lO' 7 


Table 1: Error of downstream variables at the shock between S-F and linear 
theory. Error is computed on the symmetry axis. Both absolute (a.e.) and 
relative errors (r.e.) are shown. 


( r _ o.2) the relative error between the linear theory prediction and the S-F computation 
is already 24% for entropy, suggesting strong nonlinear effects. The data also suggests that 
pressure is less affected by high e than are entropy and normal velocity. Results in Table 2 
also suggest that for e = 2 and beyond, cubic nonlinearities seem to come into play because 
the relative error is growing slightly faster than quadratically with the vortex strength. Note 
that at M - 2, flow reversal occurs at e « 3.7 which is beyond the nonlinear threshold (see 
Section 5.2). In all the cases considered, the linear results underestimate the magnitude of 
the nonlinearly computed perturbations of s 2 , p^t an< ^ u 2r - 

Further detailed checks of the accuracy of the S-F scheme have been performed in the 
two-dimensional context (Erlebacher el al. 1996). The numerical algorithm is identical to 
the current one, except for modifications required to accomodate the cylindrical geometry. 

After the successful comparison between linear theory and S-F results, we compute the 
vortex/shock interaction with ENO and S-F schemes and compare the solutions in the down- 
stream domain, specifically for the case of M — 2 and t = 2. We initialize the vortex ac- 
cording to IC-c (2.11) with q = 2, and /3 = 6. The shock is initially located at x = 0. For 
the S-F simulation, the upstream domain extends to x = 8, and the downstream domain 
extends to x = —10. with the buffer domain starting at x = —9. Parallel to the shock, 
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3 

2 

1 

0.1 

T 

0.2018 

0.1345 

0.06726 

0.006726 

entropy S 

24.05% 

9.95% 

2.37% 

0.024% 

pressure p 

4.38% 

1.90% 

0.49% 

0.005% 

velocity u n 

13.17% 

5.69% 

1.33% 

0.014% 


Table 2: Relative error of downstream variables at the shock between S-F and 
linear theory. Error is computed on the symmetry axis. 

the grid stretching parameters are A = 2.5, r max = 20 (see Equation 3.5). In computational 
space, the grid is uniform in both upstream and downstream domains, but with different grid 
spacings in the two domains. There 60 axial grid points in the upstream domain. Coarse 
(100 x 48) and fine (200 x 96) downstream grids produce almost identical contour plots of 
u ' x , indicating a converged solution. 

Conditions for the ENO scheme are as follows. The physical domain is defined by the 
region x £ [—10,8] and r £ [0,5], with a uniform grid in both directions. Coarse and fine 
resolutions grid densities are 100 x 50 and 200 x 100 respectively. Although the results from 
the two algorithms were compared at various times, we present only the results at t = 6. 
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Figure 3: Superposition of ENO versus S-F in the region x G [-4,0], r G [0, 1.5]. Contour levels of 
u x range from -9.0 to -0.6 by increments of 0.1. Solid lines: S-F, dashed lines: ENO with 200 x 100 
resolution, dotted lines: ENO with 100 X 50 resolution. 



Figure 4: Cut of Figure 3 at x = -1.9 and x = -6, r = 0.5 and r = 1.5 to show quantitative point 
by point comparison. 


Figure 3 compares contour plots of axial velocity obtained from ENO and S-F in the region 
x ^ [_5 5 _1] an d r G [0,3]. There are 4 equally spaced contour levels ranging from —9.0 to 
-0.6, with a maximum of -0.53 near x = -0.3. With a freestream velocity U\ = -0.887, the 
flow has not yet reversed. (The flow in the downstream domain is defined to have suffered 
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reversal, i.e., incipient vortex-breakdown, if the axial velocity points upstream at least at one 
point.) According to Table 3, the circulation must approach the value 4 for breakdown to 
take place. The solid line represents the S-F results which may be considered as converged 
solutions. As expected, the ENO solution approaches the S-F solution as the grid is refined. 
When the axial grid resolution is doubled, the ENO contours lie approximately half way 
between the S-F contours and the ENO contours with a coarser grid. Cuts through the 
contours of Figure 3 at r = 0 and r = 1 are shown in Figure 4 to quantify the differences 
between the solutions. The refined ENO solution is approximately half way between the 
coarse ENO and SF. This is consistent with a first order accurate solution. As explained by 
Casper & Carpenter (1995), the first order accuracy is the result of the propagation of the 
first order error near the shock through the downstream characteristics. 

5.2 Vortex breakdown regime 

The comparisons performed in the previous section demonstrate that linear theory accurately 
predicts the phenomena associated with the interaction a shock wave with an axial vortex for 
the range of vortex circulations t and shock Mach numbers M such that eM < 0(1). As the 
vortex strength is increased beyond this linear limit, the deceleration of the axial flow across 
the shock increases the likelihood of the vortex breaking down. This is best understood by 
considering the angle formed by the streamlines and the symmetry axis, (which we called in 
the Introduction as helix angle), 

0 = arctanr = arctan U6max (5.1) 

«oo 

where u ^ is the freestream axial velocity. This equation shows that 0 increases if either 
tioo decreases (stronger shock), or ue ma x increases (stronger vortex). In the case of the 
incompressible axisymmetric vortex it is known that breakdown is highly probable when 
0 reaches a critical value around 57° (Spall et al. 1987). Experiments (Cattafesta 1992; 
Delery et al. 1984), and numerical simulations (Delery et al. 1984) indicate that the vor- 
tex strength r required to induce flow reversal downstream of the shock decreases with 
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increasing Mach number. We are interested in reproducing these results with an accurate 
time-dependent simulation. 

Although both ENO and S-F are adequate to trace out the envelope of incipient break- 
down, we opted to use ENO for expediency. Spot checks at M = 2 and M = 5 with S-F 
confirm the ENO results. The initial conditions used for this were IC-b. 

The computational domain for ENO is x £ [-4,2], r £ [0,5], with a uniform grid of 
100 x 75 points.' A uniform axial vortex is placed in the upstream domain x > 0. For each 
shock Mach number, several runs are conducted with values of e on either side of its critical 
value. To this end, we compute the maximum axial velocity in the downstream domain. If 
this value is greater than zero, we consider that the flow has reversed. In the absence of 
reversal, the minimum axial velocity eventually oscillates around a constant negative value, 
wfiich we take as an indication that the solution has reached a quasi-steady state. 

The results of the parametric study are recorded in Table 3 and Figure 5. Table 3 
shows a range of Mach numbers from 1.3 to 10. The value of e for which the incipient 
breakdown occurs is given in the second column, and is the average of the two numbers in 
the third column. These are the two closest values of e we tested which are on either side 
of the breakdown curve. Clearly, determining the critical value of t to arbitrary precision 
is naturally difficult and prohibitively expensive, as is the case with physical or numerical 
experiments. Several refined ENO simulations were performed on a 200 x 150 grid to confirm 
the results obtained on the coarser mesh. 

The variation of r = r x versus Mach number is shown in Figure 5. As expected, the 
upstream vortex strength necessary to initiate breakdown decreases with increasing Mach 
number. Intuitively, this makes sense since a stronger shock leads to a stronger deceleration 
of the axial flow, thus pushing the helix angle towards its critical value. An alternative 
way to quantify the relationship between shock and vortex strength is to plot the degree of 
nonlinearity of the shock- vortex interaction versus Mach number (Fig. 5). A good curve fit 
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M 

c 

range of t 

T 

1.3 

4.93 

(4.90, 4.95) 

o.sio ; 

1.5 

4.53 

(4.50, 4.55) 

0.406 

1.7 

4.18 

(4.15, 4.20) 

0.331 

2.0 

3.78 

(3.75, 3.80) 

0.254 

2.5 

3.33 

(3.30, 3.35) 

0.179 

3.0 

3.02 

(3.00, 3.03) 

0.135 

3.5 

2.73 

(2.70, 2.75) 

0.105 

4.0 

2.64 

(2.63, 2.65) 

0.0888 j 

4.5 

2.51 

(2.50, 2.51) 

0.0750 

5.0 

2.51 

(2.50, 2.51) 

0.0675 

6.0 

2.28 

(2.25, 2.30) 

0.0511 

7.0 

2.18 

(2.15, 2.20) 

0.0419 

8.0 

2.16 

(2.15, 2.16) 

0.0363 

9.0 

2.11 

(2.10, 2.12) 

0.0315 

10.0 

2.08 

(2.05, 2.10) 

0.0280 


Table 3: The vortex breakdown or flow reversal values. 


is provided by the straight line 

eM = 1.7 M + 4.0 (5.2) 

over the range of M G [1, 10]. Above M = 6, the linear approximation tM = 1.8 M + 3.2 
provides a better fit to the data. Below this Mach number, the latter equation slightly 
underpredicts the numerical values. Equation (5.2) could serve as a practical curve for 
experimentalists and engineers. 

As explained in Delery (1984), the presence of the shock decelerates the mean upstream 
vortex thus increasing 0. A rule of thumb for incompressible flow is that when 0 reaches 
57°, the vortex is susceptible to an axisymmetric breakdown. A simple model of vortex 
breakdown was given in Cattafesta (1992) who expressed the downstream r 2 as a function 
of upstream quantities to obtain 

T2_Ut±_ p2l . . 

_ — • 2 / , 2 I W-'- 5 / 

T 1 U x2 p 2 \ sin xp + cos^ V 

where u x ( i = 1,2) are axial velocities upstream and downstream of the shock respectively, 
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Figure 5: Plot of T\ = r and tM versus Mach number for M € [1.5, 10]. Results based on ENO. 
and ip is the angle between the shock and the shock normal (Figure 1). The density ratio 


cross the shock is 


P2l — 


(7 + 1 )A/ 2 cos 2 xp 


(7 — 1 )M 2 cos 2 tp + 1 v ' 7 

Cattafesta (1992a) assumed that the shock was normal {ip = 0°) and determined a best 
constant value of 77 over the range of Mach numbers considered by Delery et al. (1984). 
This yielded t 2 = 0.6, or 0 2 ~ 30° which is at variance with the critical 0 2 = 57° for the 
incompressible vortex. Assuming that the incompressible limit for 0 2 should be reached 
at Mi leads to the conclusion that the critical value of 0 2 at which breakdown occurs 


is a function of the upstream shock Mach number. Therefore, we compute 0 2 based on 
Equation (5.3), with the values of ti taken from Table 3 to obtain the variation of 0 2 
with M shown in Figure 6. Note that the curve is approximately linear for M < 3.0. 
Extrapolated to M = 1, the 0 2 curve predicts a critical angle 0 2 = 55°, consistent with 
incompressible results. In the limit M = 1, the two curves join, and 0i = 0 2 ~ 55°. At 
higher Mach numbers, the critical angle 0 2 decreases down to 10°. Delery et al. (1984) 
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obtained experimental points on the reversal curve for M = 1.7, 2.0, 2.3. They obtained 
0 2 = 47°, 40° and 30.5°, respectively. Based on our observations relating the effects of the 
pressure core, and the resulting shock motion, the discrepencies could be partially explained 
by the presence of an axial velocity deficit present in the experiments, but absent from the 
simulations. 



Figure 6: 0i and ©2 versus Mach number. 

If t is increased beyond the point where the vortex breaks down, the shock eventually 
forms a triple point, with a more complicated structure in the recirculating region. This 
regime cannot be computed by the S-F code because of the appearance of internal shocks in 
the downstream This regime is now studied exclusively with the ENO algorithm. 

5.3 Shock bifurcation regime 

To complete the study of flow configurations as a function of vortex strength and shock 
Mach number, we perform ENO calculations for e = 7 and 9, at M = 2 and 4. The physical 
domain extends from x = — 4 to x = 8 in the axial direction, and to r = 5 radially. Grid 
resolution is 450 x 150. These parameters were chosen to ensure the occurrence of a triple 
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point in the shock, hereafter referred to as shock bifurcation. The triple point is expected to 
form at r = 0(1) since the vortex core radius is the only physical length scale in the problem. 
The chosen values of c far exceed the breakdown criterion; the flow is highly nonlinear. 

We consider two different set of initial conditions — IC-a and IC-c. Recall that IC-a 
refers to a vorticity profile and perturbation temperature constant across the upstream and 
downstream domains; whereas IC-c has zero vorticity downstream and an upstream vorticity 
multiplied by a spatial blending function (Section 2). The structure of the flow is a function 
of the initial conditions (compare figures 8 and 11). 

Contour plots of |Vp| and |Vp| for a M = 2, c = 7 flow with IC-c at t = 11 are shown in 
Figure 7. For comparison, density contours are shown in Figure 8. The gradient functions 
are normalized so that their maximum is unity. Contours in the range [0.03,0.04] are plotted. 
The choice of variables is consistent with the features visible in a typical Schlieren photo. The 
flow slip lines are easily identified from the features present in the density gradient contours, 
and absent from the pressure gradient plot (Pi, P 2 , and H in Figure 7). The upstream 
pressure is minimum in the vortex core, specifically at the axis, and increases radially to 
its freestream value. As the plane shock interacts with the vortex (at Pi), the local shock 
Mach number is therefore higher in the core region than elsewhere. In other words, the 
shock bulges forward in the vicinity of the axis. As the obliqueness of the primary shock (A) 
increases, a local region of supersonic flow eventually forms downstream of (A), necessarily 
terminated by another shock ( B ) to match with the downstream subsonic flow, thus forming 
a lambda shock with the so-called triple point (P 2 ). Another weak shock (C) is evident in the 
supersonic region, between shocks (A) and ( B ). The sliplines ( P t ) and (E 2 ) are separated by 
the foot of the shock (C), and they coincide with the Mach=l contour line which emanates 
from the point of intersection of the primary shock (A) with the axis of symmetry (Pi). 
Thus, the supersonic region is delimited by the shocks (A), ( B ) and the sonic lines (Pi) and 
(p 2 ). As the primary shock (A) begins to move upstream in the region near the axis, an 
azimuthal vortex ring is formed. This vortex ring, which has an elliptical cross section in the 
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Figure 7: Schematic of flow configuration for a smooth start case. 
x — r plane whose major axis slowly rotates counterclockwise, moves steadily away from the 
axis. The flow in the neighborhood of the elliptic vortex closest to the center is supersonic, 
while the rest is subsonic. As the flow turns supersonically within the vortex (Figure 9), it 
is decelerated through the shocks ( B ) and the triple point structure at the tip of the normal 
shock ( D ). The velocity vectors near the axis in Figure 9 suggest a structure analogous to 
that of an expanding nozzle with a shock consistent with subsonic exit conditions. 
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- 5 - 3-1 1 3 5 


x 

Figure 9: Velocity vector plots in the x — r plane, superimposed on contours of Mach number 
(projected in the x — r plane). Contours are at intervals of 0.5. The vectors are plotted every 10 
points in x and every 4 points in y . 



Figure 10: Position of P\ and P 2 as a function of time. The initial data and grid is given in the 
text. 
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Figure 8 shows density contour plots at t = 3.7,11 for the IC-c, with a spatial ramp 
function defined by Equation (2.11) (a = 2, (3 = 6). Although we do not present the results 
here, we have also simulated the abrupt start case defined by IC-b, and find that the results 
are qualitatively similar to Figure 8, but with a time shift. Figure 8a ( t = 3) shows the 
density prior to the formation of the triple point. The triple point forms around t = 4. 
By t = 7, the triple-point shock configuration is well formed (Figure 8b), and as evidenced 
by Figure 8c at t = 11, this formation seems to evolve self-similarly. The structure of the 
flow in the recirculatig region is presented in Figure 9 in which vector plots of the flow in 
the x — r plane are superimposed on contours of Mach number. The contours range from 
0 to 2.5 in increments of 0.5. The results are shown at t = 11. As the flow begins to wrap 
around the subsonic region (E), it decelerates, and becomes subsonic. At this point, the 
flow reaccelerates, eventually becomes supersonic, and is eventually slowed down by passage 
through a normal M = 2 shock (D). 

To test whether the flow maintains its self-similar motion for longer time periods, we first 
coarsened the grid to 100 x 50 on the same physical domain and compare the primary shock 
structure to the finer grid results of Figure 8. We find, as expected, that the coarseness of 
the grid does not modify the characteristics of the downstream flow. Therefore, we increased 
the size of the physical domain to x G [-9, 15], r G [0, 12], and ran the simulation to t = 28. 
At this late stage, the structure of the primary and secondary shocks are still the same as 
for t — 11, but spatially enlarged. We track the position of Pi (the point of intersection 
of the shock and the axis) and P 2 (the shock triple-point) and find that P 2 is moving at 
a constant velocity towards the freestream (Figure 10). On the other hand, the speed of 
Pi seems to slow down as it moves upstream. These results suggest the possibility of a 
theoretical framework within which a self-similar flow emerges for particular combinations 
of the parameters M and t. 

In contrast to results with the initial conditions IC-c, Figures lla-c show the time evo- 
lution of the density for case IC-a, where a vortex is superimposed initially on a shocked 


30 






uniform flow. Clearly, the evolution appears to be self-similar, somewhat like that in Fig. 8, 
but evolving on a slower time scale. Rough velocity measurements of points Pi and P 2 be- 
tween t = 7 and t = 11 indicate that for IC-a, Pi moves approximately 3.5 times slower, and 
P 2 moves about 2.5 times slower, than the corresponding points of case IC-c. This indicates 
that the shock is more oblique when the vortex interacts with the shock abruptly. A physical 
explanation is offered by considering the initial pressure jump across the shock, which for 
our choice of parameters, is about twice as strong in case IC-b. This indicates that case IC-b 
generates a shock with stronger acceleration into the upstream domain near the symmetry 
axis. 

Once the triple point is formed and starts to move towards the freestream, what is its 
ultimate destiny? Does it continue to move at a constant velocity, or does it slow down, 
and eventually stop, setting the stage for a possible steady state solution? We try to answer 
this question by appealing to dimensional analysis. The physical parameters that govern 
this problem are the mean upstream velocity Ui, the upstream mean sound speed Ci, the 
normalized vortex strength t = T y/^j (r 0 c*). the vortex core radius r 0 , the spatial coordinates 
x* and c*, and time t*. (An asterisk denotes dimensional quantities.) The downstream 
velocity field therefore takes the general form 


n , , ... _ 3 -* V * Cit m 

r 0 r 0 r 0 


(5.1) 


In the limit of vanishing core radius, or equivalently, as — > oo, the velocity field must 

take the self-similar form 


“K^r.TT )=c7(M.<. 1 


y 


-J- 


(5.2) 


f t*' ' ' ' ' Cit * ’ c t f 

Expressions for the velocity of the triple point P 2 and that of the intersection of the 
primary shock and the axis ( P \ ) are easily obtained from (5.1 ) by suppressing the dependence 
on the spatial coordinates: 

(5.3) 


ci /(M,e, — ) 


r o 
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which implies a constant velocity as the third argument becomes very large. Of course, this 
velocity might be zero. 

Numerical results seem to indicate that there is a region of (e, M) parameter space for 
which the function / is independent of time. However that need not be the case. A theoretical 
study of this matter would prove most interesting. For example, what is the structure of / 
as t — >• oo, T — > oo, or M -4 co? Intuitively, one could argue that once out of range of the 
upstream vortex filament, that having reached a constant velocity, the triple point would 
continue on forever. Given the assumption that the point Pi moves at constant velocity, a 
slowdown of the triple point would increase the obliqueness of the shock, eventually leading 
to its breakdown. Simulations at higher values of t have led to shock breakdown. 

To understand better the effect of M and e on the triple point motion, we perform and 
display results at M = 2 and 4, for e = 7 and 9. These are presented in Figure 12. Comparing 
Figures 12a (M = 2,e = 7) and 12b (M = 2, e = 9), one sees that the velocity of the triple 
point P 2 (Figure 7) does not depend on the vortex strength. However, Pi moves at a much 
faster rate as e is increased. This is consistent with a lower vortex core pressure, and thus a 
greater shock velocity near the axis. These results also hold for M = 4 at the same values 
of e (Figures 12c-d). On the other hand, the velocity of both Pi and P 2 increases with Mach 
number. 

These results allow us to draw some conclusions relating to the form of Equation (5.3). 
The lack of dependence of vp 2 on e implies that 

vp 2 = ci f(M, — ) (5.4) 

r o 

which is only a function of Mach number when the velocity is independent of time (i.e., 
r 0 — > 0). On the other hand, vp l retains the general form (5.3). 

Now consider the limit as M — > oo for which either cj -4 0 or «i oo. We consider 
both cases. But first, consider the radial temperature profile. If one assumes that the entire 
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vortex is isentropic, Equation (2.6) implies that 


T(r) = Ti - o 2 (r)e 2 


(5.5) 


where a(r) is a positive function and cv(oo) = 0. Therefore, the Mach number in the vortex 
core is higher than in the freestream. Its radial profile is 



which gives a ratio of core to freestream Mach number 


Mo _ 1 

Mi ~ (1 -a 2 (0 )7 c 2 /c 2 ) 1 / 2 ' 


(5.7) 


If U\ -4 oo, = /(e,Cj) with t < to insure a positive core temperature. Therefore, 

the Mach number and vortex strength are independent parameters. If, on the other hand, 
Ci -4 0, the ratio Mq/Mi is still given by (5.7), except that the maximum vortex strength 
decreases proportionally to Cj. In this case, one cannot take the limit Mi -4 oo while 
keeping e fixed. In other words, the aforementionned expressions for the triple points , and 
their “self-similar’ properties, are found for increasingly weaker vortices as the Mach number 


is increased. 


Combining the above statements, as M\ -4 oo, (ci -4 0), the velocity of Pi and P 2 satisfy 
either vp t = U\ f(c,Uit*/r 0 ) or vp t — Ci f{^,Uit*/r 0 ). In the latter case, the velocity drops 
to zero in the limit M\ -4 00 . We assume that e is kept fixed as ci — > 0, e.g., by decreasing 
the vortex circulation. Since t < O(ci), one must have T < 0(c 2 ). 

When Mi — » 00 (Ui -4 00), vp t = C\ /(c, Ci<*/r 0 ), which implies that the triple point has 
a velocity vp 2 = Ci /(e, Cif*/r 0 ). In the above discussion, we have assumed that the generic 
function /(•••) remained finite as its various limits were taken. However, this need not be 
true. If this assumption fails, more general forms must be considered. 

Figures 13a-13d show the solution of the IC-c case at t = 11 for p. u x , ue , and axial 
vorticity w x ( M = 2, e = 7). An “airfoil”-like structure is apparent in Figure 13b. This 
region acts as a solid body around which the flow rapidly accelerates (Figure 13c). If the 
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acceleration is sufficiently strong, the flow becomes supersonic as it redirects itself in the 
upstream direction. To adjust back to the subsonic field adjacent to the primary shock, a 
secondary shock is formed at the axis (D in Fig. 7). The secondary shock forms an almost 
perfect normal Mach disk. The axial location of this shock is approximately steady and 
has a Mach number of 1.34. Note that this recirculating zone does not appear in pressure, 
but is present in density plots, which indicates a region of localized entropy variation. As 
expected, the axial vorticity is continuous across the shock, but has strong variations inside 
the “airfoil”. This is consistent with strong variations of vg in this region. 

6 Concluding Remarks 

In this paper, we have studied the interaction of a shock with a longitudinal isentropic 
vortex over a wide range of Mach numbers and vortex strengths. Three regimes are brought 
to light. In the first, the vortex strength t satisfies eM < 1, in which case, linear results are 
valid. As e increases, the nonlinear effects increase. First, a flow reversal occurs downstream 
of the shock, accompanied by a vortex breakdown. This effect is primarily due to the 
deceleration of the mean flow across the shock, which decreases the helicity of the vortex, 
thus eventually leading to breakdown. We determined numerically e as a function of Mach 
number for which incipient vorticity breakdown occurs. When the ordinate is expressed as 
cM , this curve becomes approximately linear, particularly at the higher Mach numbers. As t 
is increased even further, the resulting obliqueness of the shock leads to a supersonic pocket 
in the downstream region which readjusts to subsonic conditions through a secondary shock 
emanating from the primary shock at the so-called triple point. A parameter study in this 
range of vortex strengths using ENO, uncovered a regime in wffiich the motion of the triple 
point and of the point of intersection of the primary shock with the axis of symmetry looks 
self-similar. Further study uncovered a regime in which the triple point velocity, normalized 
by Ci, was only a function of the shock Mach number. Under extreme conditions, the shock 
can even break, when a supersonic upstream flow can no longer be maintained. This is due 
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to the strong Mach number gradient between the vortex core and its edge. This condition is 
easiest to reach for relatively weak shocks. Currently, there are no nonlinear theories which 
support these results. Theoretical developments, particularly those based on self-similar 
notions, might go a long way towards understanding in more detail some of the mechanisms 
described herein. 
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